%% 多元分析Multivariate analysis
%% 1.聚类分析Cluster analysis
% 聚类分析又称为群分析，是对多个样本或指标进行定量分类的一种多元统计分析方法
% 。对样本进行分类称为Q型聚类分析，对指标进行分类称为R型聚类分析。下面通过实
% 例来介绍两类聚类分析。
clear,clc
load data.mat % 该数据为某地区6项空气质量指标在3126天的变化情况
n=300;
shuju1=shuju(1:n,:);% 取前300项样本进行Q型聚类分析
shuju1=zscore(shuju1);% 数据标准化
y=pdist(shuju1,"euclidean");% 计算各行向量两两之间距离，常见方法如下：
% 'euclidean'欧氏距离(默认)；'seuclidean'标准欧氏距离；'cityblock'绝对值距离；
% 'minkowski'闵氏距离，后边跟数字k表示q=k时距离；'chebychev'切比雪夫距离；
% 'mahalanobis'马氏距离；'cosine'夹角余弦；'correlation'样本相关系数；
% 'spearman'spearman秩相关系数。
% 若行数为k，其返回为1列向量，元素个数为组合数C(k)(2)。
y1=squareform(y);% 将y转化为矩阵形式，矩阵中元素(i,j)对应原始数据集对象i和j
% 间距离。
z=linkage(y,"average");% 生成聚类树，常见方法如下：
% 'single'最短距离(默认)；'average'无权平均距离；'centroid'重心距离；
% 'complete'最大距离；'median'赋权重心距离；'ward'离差平方和方法；
% 'weighted'赋权平均距离。
% 输出为包含聚类树信息的(m-1)*3矩阵。
h=dendrogram(z);% 画聚类图
t=cluster(z,"maxclust",5);% 创建聚类，分为k类
% R型聚类只需将原始数据进行转置后分析，基本分析方法与Q型聚类相同，此处不再赘述。
%% 2.主成分分析Principal component analysis
% 主成分分析可以将相关性高的变量转化为彼此相互独立的变量，是一种数据降维方法。
% 在变量进行最小二乘估计时，当变量间相关性较高时，可能表现出不稳定性。此时可
% 利用主成分分析方法将原来自变量变换为另一组主成分变量，用新的变量进行回归分
% 析。
clear,clc
load data.mat % 该数据为某地区6项空气质量指标在3126天的变化情况
shuju1=shuju';% 分析六项空气指标间主成分
shuju1=shuju1(:,1:300);% 取300项样本进行分析
shuju1=zscore(shuju1);% 数据标准化
y1=corrcoef(shuju1');% 计算相关系数
[vec1,lamda,rate]=pcacov(y1);% 主成分分析，lamda为相关系数矩阵的特征值，rate
% 为各个主成分的贡献率
contr=cumsum(rate);% 计算累积贡献率
f=repmat(sign(sum(vec1)),size(vec1,1),1);% 构造与vec1同维数，元素为+-1的矩阵
vec2=vec1.*f;% 修改特征向量正负号，使得每个特征向量分量和为正
num=4;% 选取主成分个数
df=shuju1'*vec2(:,1:num);% 计算各个主成分得分
tf=df*rate(1:num)/100;% 以num个主成分贡献率作为权重，构建综合评价模型，可计
% 算各样本的综合得分。
%% 3.因子分析Factor analysis
% 因子分析可以看成主成分分析的推广，它也是多元统计分析中常用的一种降维方式，
% 因子分析把方差划归为不同的起因因子，需要构造因子模型，通过潜在的假象变量和
% 随机影响变量的线性组合表示原始变量。
% 因子分析主要有以下步骤：
% (1)选择分析的变量。因子分析的前提条件时观测变量间有较强的相关性，因为如果
% 变量之间无相关性或相关性较小，他们之间就不会有共享因子。
% (2)计算所选原始变量的相关系数矩阵。
% (3)提出公共因子。这一步需要确定因子求解的方法和因子的个数，按照因子累积方差
% 贡献率一般要达到60%才能符合要求。
% (4)有时提出的因子很难解释，需要通过坐标变换使每个原始变量在尽可能少的因子之
% 间有密切的联系，这样因子解的实际含义更容易解释。
% (5)计算因子得分。
% 本节通过实例说明因子分析的过程：
clear,clc
load data.mat % 该数据为某地区6项空气质量指标在3126天的变化情况
shuju1=shuju';% 分析六项空气指标间主成分
shuju1=shuju1(:,1:300);% 取300项样本进行分析
shuju1=zscore(shuju1);% 数据标准化
y1=corrcoef(shuju1');% 计算相关系数
% 估计因子载荷矩阵：
% (1)主成分法估计
[vec1,lamda,rate]=pcacov(y1);% 主成分法
f1=repmat(sign(sum(vec1)),size(vec1,1),1);% 构造与vec1同维数，元素为+-1的矩阵
vec2=vec1.*f1;% 修改特征向量正负号，使得每个特征向量分量和为正
f2=repmat(sqrt(lamda)',size(vec2,1),1);% 构造与vec2同维数的矩阵
a=vec2.*f2;% 构造全部因子的载荷矩阵
num=3;% 提取num个主因子的载荷矩阵
aa=a(:,1:num);
tcha2=diag(y1-aa*aa');% 计算num个因子的特殊方差
ccha2=y1-aa*aa'-diag(tcha2);% 计算num个因子时残差矩阵
% (2)主因子法估计
nk=size(y1,1);rt=abs(y1);% 求相关系数阵所有元素绝对值
rt(1:nk+1:nk^2)=0;% 对角线元素赋值0
rstar=y1;rstar(1:nk+1:nk^2)=max(transpose(rt));
% 将y1换为rstar，该步之后进行与主成分法估计相同的步骤。
% (3)最大似然估计
% 当y1正定时，可以使用该方法[lambda,psi]=factoran(y1,1,'xtype','cov');
% 因子旋转：
% 因子旋转的目的是使因子载荷矩阵结构简化，使载荷矩阵每列或行的元素平方值向0
% 和1两级分化。主要有三种正交旋转法：方差最大法、四次方最大法和等量最大法。
[lambda2,t]=rotatefactors(aa,"Method","varimax");
%% 4.判别分析Discriminant analysis
% 判别分析是根据所研究的个体的观测指标来推断该个体所属类别的一种统计方法，在
% 自然科学和社会科学的研究中经常会遇到这类问题。判别分析问题用统计语言来描述，
% 就是已知有q个总体X1,X2,...,Xq，它们的分布函数分别来自F1(x),F2(x),...,Fq(x)
% 每个分布函数都是p维函数，对于给定的样本x，要判断它来自哪一个总体。
% 基本的判别方法有距离判别、Bayes判别，Fisher判别和距离判别。在机器学习章节对
% 于分类问题的研究方法进行了细致讨论，本节介绍距离判别和Fisher判别方法。
% 距离判别主要结合已知样本的均值和协方差信息定义距离判别函数，判断待测点距离
% 哪类样本更近。
% 判别函数：w(x)=(x-u2)'/COV2*(x-u2)-(x-u1)'/COV1*(x-u1);
% 式中：COVi=1/(ni-1)*sum((xi-ui).^2);ui为第i类数据的均值，ni为第i类数据的个
% 数。则当w(x)>=0时，待测点属于1类，否则属于2类。
% Fisher判别的基本思想是投影，即将表面上不易分类的数据通过投影到某个方向上，
% 使得投影类与类之间得以分离的一种判别方法。Fisher判别要求有公共的协方差矩阵
% COV，则有Fisher判别函数：
% w(x)=[x-1/2*(u1+u2)]'/COV*(u1-u2);当w(x)>=0时属于1类，否则属于2类。
clear,clc
% 设有原始数据：
n=1200;% 观测次数
m=4;% 自变量个数
x=-20+40*rand(n,m);% 观测数据为m项，每项n个数据，进行以下分析时，一般假定各项数
% 据间相关性较小。
y=zeros(n,1);
w=[-1,-1,1,1];b=0;% 为了让数据尽可能线性可分，先预设w,b设计原始数据
for i=1:n
    y0=randn(1)+sum(w.*x(i,:))+b;% 利用高斯噪声和线性值设计y1，理论上beta
    if y0<0
        y(i)=1;% 对数据预分类
    else
        y(i)=2;
    end
end
% 利用matlab内置方法进行训练分类：
[x1,y1]=classify(x,x,y,"mahalanobis");% 使用距离判别，返回x1为预测值，y1为误
% 判率。输入第一个x为待测样本，第二个x为已知样本，y为已知样本响应。分类方法主
% 要有：'linear'线性；'quadratic'二次；'mahalanobis'马氏距离
%% 5.典型相关分析Canonical correlation analysis
% 通常情况下，为了研究两组变量[x1,x2,...,xp],[y1,y2,...,yq]的相关关系，可以
% 分别计算两组变量之间的全部相关系数，共有p*q个相关系数。典型相关分析分别找出
% 两组变量各自的某个线性组合，讨论线性组合之间的相关关系。
clear,clc
load data.mat % 该数据为某地区6项空气质量指标在3126天的变化情况
shuju1=shuju';% 分析六项空气指标间主成分
shuju1=shuju1(:,1:300);% 取300项样本进行分析
shuju1=zscore(shuju1);% 数据标准化
y1=corrcoef(shuju1');% 计算相关系数
% 现在利用典型相关分析6项空气指标中前3项与后3项的相关关系：
n1=3;n2=3;num=min(n1,n2);
s1=y1(1:n1,1:n1);s12=y1(1:n1,n1+1:end);
s21=s12';s2=y1(n1+1:end,n1+1:end);
m1=inv(s1)*s12*inv(s2)*s21;m2=inv(s2)*s21*inv(s1)*s12;
[vec1,val1]=eig(m1);% 求m1特征向量和特征值
for i=1:n1
    vec1(:,i)=vec1(:,i)/sqrt(vec1(:,i)'*s1*vec1(:,i));
    % 特征向量归一化
    vec1(:,i)=vec1(:,i)*sign(sum(vec1(:,i)));
end
val1=sqrt(diag(val1));% 计算特征值的平方根
[val1,ind1]=sort(val1,'descend');% 从大到小排列特征值
a=vec1(:,ind1(1:num));% 第一组系数阵
dcoef1=val1(1:num);% 提出典型相关系数
[vec2,val2]=eig(m2);% 求m1特征向量和特征值
for i=1:n2
    vec2(:,i)=vec2(:,i)/sqrt(vec2(:,i)'*s2*vec2(:,i));
    % 特征向量归一化
    vec2(:,i)=vec2(:,i)*sign(sum(vec2(:,i)));
end
val2=sqrt(diag(val2));% 计算特征值的平方根
[val2,ind2]=sort(val2,'descend');% 从大到小排列特征值
b=vec2(:,ind2(1:num));% 第一组系数阵
dcoef2=val2(1:num);% 提出典型相关系数
%% 6.对应分析Correspondence analysis
% 当研究变量相关关系时采用R型因子分析，研究样品间相互关系时则采用Q型因子分析
% 但无论是R型还是Q型因子分析都未很好的揭示变量和样品间的双重关系。对应分析在
% 因子分析的基础上发展而来，它对原始数据采用适当的标度方法，把R型和Q型分析结
% 合起来，同时得到两方面结果———在同一因子平面上对变量和样品仪器进行分类，从而
% 解释所研究样品和变量间的内在关系。
% 本节通过示例说明对应分析过程。
clear,clc
load data.mat % 该数据为某地区6项空气质量指标在3126天的变化情况
shuju1=shuju;% 分析六项空气指标间主成分
shuju1=shuju1(1:300,:);% 取300项样本进行分析
[num1,num0]=size(shuju1);
for i=1:num0
    shuju1(:,i)=(shuju1(:,i)-min(shuju1(:,i)))/(max(shuju1(:,i))-min(shuju1(:,i)));
    % 数据归一化至(0,1)
end
% (1)化数据阵为规格化"概率"矩阵P：
T=sum(sum(shuju1));
P=shuju1/T;% 求概率矩阵P
r=sum(P,2);c=sum(P)';% 求边缘概率
% (2)进行数据对应变换：
B=zeros(num1,num0);
for i=1:num0
    for j=1:num1
        B(j,i)=(P(j,i)-r(j)*c(i))/(r(j)*c(i));
    end
end
% (3)计算行轮廓分布：
R=zeros(num1,num0);
for i=1:num1
    R(i,:)=P(i,:)/r(i);
end
Dr=diag(r);Dc=diag(c);
% (4)计算列轮廓分布：
C=zeros(num1,num0);
for i=1:num0
    C(:,i)=P(:,i)/c(i);
end
% (5)对B做奇异值分解：
[u,s,v]=svd(B,'econ');
w=sign(sum(v));vb=v.*w;ub=u.*w;
lamda=diag(s).^2;% 计算B'*B特征值，即计算主惯量
ksi=T*(lamda);% 计算卡方统计量的分解
T_ksi=sum(ksi);% 计算总卡方统计量
con_rate=lamda/sum(lamda);% 计算统计量
cum_rate=cumsum(con_rate);% 计算累积贡献率
beta=diag(r.^(-1/2))*ub;% 求加权特征向量
G=beta*s;% 求行轮廓坐标
alpha=diag(c.^(-1/2))*vb;
F=alpha*s;% 求列轮廓坐标F
%% 7.多维标度法Multidimensional scaling method
% 在实际中往往会碰到这样的问题：有n个有多个指标(变量)反映的客体，但反映客体的
% 的指标个数是多少不清楚，甚至指标本身是什么也是模糊的，所能知道的仅仅是这n个
% 客体之间的某种距离或者某种相似性，我们希望仅由这种距离或相似性出发，在较低
% 维欧氏空间中把n个客体图形描绘出来，从而尽可能揭示这n个客体之间的真实结构关
% 系。
% 本节通过实例给出求多维标度法经典解的步骤。
clear,clc
D=[0,1,sqrt(3),2,sqrt(3),1,1;0,0,1,sqrt(3),2,sqrt(3),1;0,0,0,1,sqrt(3),2,1
    0,0,0,0,1,sqrt(3),1;0,0,0,0,0,1,1;0,0,0,0,0,0,1;0,0,0,0,0,0,0];
d=D+D';% 构造距离矩阵
[y,eigvals]=cmdscale(d);% 求经典解，d为实对称矩阵或pdist函数的行向量输出
plot(y(:,1),y(:,2),'o','Color','k','LineWidth',1.3)% 绘制点坐标
% 通过特征值求经典解：
A=-d.^2/2;% 构造A矩阵
n=size(A,1);H=eye(n)-ones(n)/n;% 构造H矩阵
B=H*A*H;
[vec1,val1]=eig(B);% 求B特征值和特征向量
[val2,ind]=sort(diag(val1),'descend');% 将特征按从大到小排列
vec2=vec1(:,ind);% 将特征值也按以上顺序排列
vec3=orth(vec2(:,[1,2]));% 构造正交特征向量
point=vec3.*sqrt(val2(1:2)');% 利用矩阵广播求点的坐标
